Influence of elevation on the species–area relationship
Abstract
Aim
Species–area relationships (SARs) are among the best investigated patterns in ecology, yet the shape of the function that should describe SARs and the biological meaning of the function parameters are disputed. Elevational gradients offer the opportunity of investigating how biodiversity responds to large variations in environmental characteristics within small geographical areas. We asked which function describes SARs at different elevations and explored how variations in environmental characteristics influence SAR shape.
Location
Alborz Mountains (Iran).
Taxon
Vascular plants.
Methods
We used sets of nested plots (0.001 to 100 m2) placed at 100 m intervals from 2,000 to 4,500 m elevation to construct series of nested SARs as species accumulation curves. Then, we used these curves to assess the appropriateness of different SAR functions at different elevations. We investigated how parameters of the power function varied along the elevational gradient in response to variation in environmental parameters (ruggedness, temperature, precipitation, exposed rock, percentages of soil sand and total nitrogen, and productivity, expressed by the normalized difference vegetation index).
Results
The most frequently observed best fit model was the power function, which is controlled by two parameters: z (the velocity in species accumulation with sampled area) and c (the species richness per unit area). z was positively influenced by temperature and soil nitrogen, decreasing with elevation. c was positively influenced by temperature and soil nitrogen, and negatively by rock cover, decreasing with elevation.
Main conclusions
The decrease in c‐values with elevation is consistent with the altitudinal decrease in species richness and is explained by the increase in bare rock. By contrast, c was positively influenced by temperature and total nitrogen, which are two factors promoting plant growth. Similarly, z‐values decreased with elevation, thus indicating a decrease in beta diversity.
1 INTRODUCTION
The species–area relationship (SAR), that is, the increase in number of species with area, is one of the most widespread and investigated patterns in ecology (Lomolino, 2000, 2001). In fact, there are two main types of SARs: those arising from true isolates (ISARs, island species–area relationships, in which the sampled areas are isolates and the function fitted is based on how many species are found in each sampled area) and those arising from the progressive aggregation of smaller sampling areas into larger areas (species accumulation curves, SACs, which present cumulative counts of increased species number with sampling area; see Matthews, Guilhaumon, Triantis, Borregaard, & Whittaker, 2016 for the distinction between these two concepts). Scheiner (2003) distinguished several types of SARs, according to the underlying sampling schemes: (a) strictly nested quadrats (Type I curves); (b) quadrats arrayed in a contiguous grid (Type II curves); (c) quadrats arrayed in a regular but non‐contiguous grid (Type III curves); or (d) areas of varying size, often islands (Type IV curves). While Type IV corresponds to the ISARs, the other types are forms of SACs. After decades of comparative research, it appears that no single mathematical model can adequately describe the SAR in all cases, the best function for any given study system being identifiable only empirically (Connor & McCoy, 1979; Tjørve, 2009; Triantis, Guilhaumon, & Whittaker, 2012; Williams, Lamont, & Henstridge, 2009). However, the power function S = cAz, where S is the number of species, A is the area and c and z are fitted parameters, has been found to be the model that best describes the SAR in most cases, which led to its widespread use in most studies (Dengler et al., 2019; Matthews, Guilhaumon, et al., 2016; Matthews, Rigal, Triantis, & Whittaker, 2019; Triantis et al., 2012).
The biological meaning of the parameters of the various functions proposed to model the SAR is debated, with the most complex functions being those more difficult to interpret (Dengler, 2009b; Triantis et al., 2012). The meaning of the two parameters of the power function has attracted the interest of biogeographers and ecologists (Connor & McCoy, 1979; Matthews, Rigal, et al., 2019; Triantis et al., 2012). In general, c‐values express the number of species per unit area, whereas z‐values indicate how fast species accumulate with increasing area (Fattorini, Borges, Dapporto, & Strona, 2017). The power function is usually applied in its linearized form: 
, because, in this way, the function can be easily fitted using a linear regression approach and the two parameters are more easily interpreted: c is the intercept of the regression equation (i.e. the expected number of species per unit area), and z is the slope of the regression equation (i.e. how fast species increase with increasing area). It is important to note that z is not the slope for the true species–area relationship (i.e. in arithmetic space), which is influenced by both c and z (Lomolino, 2001). There is a debate around the relationship between the power function parameters and different forms of diversity. Since c is the expected number of species per unit area it can be considered a measure of α‐diversity. The possible relationship between the z‐parameter and β‐diversity is disputed, partly as a consequence of the multitude of different definitions of β‐diversity (e.g. Jurasinski, Retzer, & Beierkuhnlein, 2009; Tuomisto, 2010). Although several authors considered the z‐parameter as a measure of β‐diversity (e.g. Drakare, Lennon, & Hillebrand, 2006; Lennon, Koleff, Greenwood, & Gaston, 2001; Polyakova et al., 2016; Ricotta, Carranza, & Avena, 2002; Scheiner, 2003, 2004; Sreekar et al., 2018), the z‐value of a species accumulation curve does not quantify compositional heterogeneity, but the rate of change in α‐diversity with increasing sampling unit size (Tuomisto, 2010). Therefore, the rate of rise of a power curve is a function of both c and z, and z alone does not express β‐diversity; however, larger values of z correspond to faster accumulations of species as area increases, thus indicating higher rates of species turnover, although limited to the addition of new species.
Most SAR parameter research has investigated causes of variation in c‐ and z‐values across different spatial and temporal scales, environmental conditions and taxa (Fattorini et al., 2017). However, studies exploring systematic variations in c‐ and z‐values are rare, because they imply the construction of a number of SARs along some ecogeographical gradient. Latitude and elevation are the two most obvious ecogeographical gradients, because a number of ecological variables (especially climatic) are known to vary with them. However, few studies have investigated how SARs vary along these gradients. In fact, Connor and McCoy (1979) and Triantis et al. (2012) failed to find a latitudinal effect on z‐ or c‐values of island SARs. However, Drakare et al. (2006), using a dataset comprising SARs for a variety of organisms, habitats and locations, found that z significantly decreased with increasing latitude. Qian, Fridley, and Palmer (2007) also found that z consistently decreased with increasing latitude for the North American flora, whereas Lyons and Willig (2002), using bats and marsupials of North and South America, found that z increased with latitude, while c exhibited a reverse pattern. Thus, despite the great interest of ecologists and biogeographers in latitudinal gradients (Fattorini & Ulrich, 2012; Lomolino, Riddle, Whittaker, & Brown, 2010; Preisser, 2019), it is not clear if and how z and c vary with latitude.
Patterns of variation in biodiversity along elevational gradients have become increasingly popular in recent years because they offer the intriguing opportunity of investigating how biodiversity responds to large variations in many environmental characteristics (such as climate, land use, soil, etc.) within small geographical areas (Fattorini, Di Biase, & Chiarucci, 2019). Quite surprisingly, elevational variation in SARs has been substantially unexplored so far, the only exception being a study by Baumann, Weiser, Chiarucci, Jentsch, and Dengler (2016) of alpine grassland vegetation, in which z‐values of nested plots (from 0.0001 to 100 m2) were found to be positively influenced by elevation.
In this paper, we used plant species distribution data collected with a nested plot design along two elevational gradients in the Alborz Mountains (Iran) to investigate how SARs vary with elevation and associated environmental characteristics. To the best of our knowledge, this is the first study that explicitly relates SAR shapes with elevation and environmental variables, thus linking two of the most important biogeographical patterns of biodiversity: the SAR and the elevational gradient. Also, by analysing SAR variation with environmental characteristics, we can shed further light on the biological meaning of SAR parameters. Specifically, the aim of the present study was to answer the following questions: (a) Given the multiplicity of the functions proposed to model the SAR, which function best describes the SAR at different elevations? (b) If the same function is applied at different elevations, how do model parameters change along the gradient? (c) Which environmental variables best explain this variation in the function parameters along the gradient?
2 MATERIALS AND METHODS
2.1 Study area
The study was conducted in the central Alborz Mountains. This mountain system is located in the north of Iran, at the southern shore of the Caspian Sea, and ranges from 26 m below sea level to 5,671 m above sea level, at the summit of Mt. Damavand. In response to climatic variation along elevational gradients, mountain vegetation tends to form distinct vegetational belts at different elevations (Fattorini et al., 2019). From the Caspian Sea level up to an elevation of ca. 2,400 m, a temperate climate allows the presence of Hyrcanian forests, a form of temperate deciduous broad‐leaved forest dominated mainly by the Oriental Beech (Fagus orientalis). These broad‐leaved forests date back 25 to 50 million years and are recognized as a distinct ecoregion (Olson & Dinerstein, 1998; Olson et al., 2001). The upper limit of these forests is dominated by Quercus macranthera, which marks the tree limit (Noroozi & Körner, 2018). Over the tree limit (which is around 2,400–2,800 m), there is a transition towards the subalpine belt, which is characterized by a colder and subhumid climate, and which is followed by an even colder and dryer alpine belt. The subalpine and alpine belts, which typically extend between 3,000 and 3,900 m (Klein, 1994), are characterized by the presence of grasslands mixed with thorny cushion shrubs, meadows and mountain pastures. The subnival–nival belt (ca. 3,900–4,800 m) is mainly covered by scree plant formations (Frey & Probst, 1986; Klein, 1982; Noroozi, Willner, Pauli, & Grabherr, 2014).
2.2 Data sampling
The study was focused on grassland vegetation. Sampling was done from the upper edge of the Hyrcanian forest belt to the subnival–nival belt along two elevational transects: the first transect (transect R) was placed on the Mt. Rostam‐Nisht (between N 36° 26′ 16.1″, E 051° 03′ 23.2″ and N 36° 24′ 05.9″, E 050° 57′ 43.7″) from 2,000 to 4,500 m, and the second (transect L) on Mt. Lashgarak the Great (hereafter Mt. Lashgarak, between N 36° 19′ 21.2″, E 051° 02′ 31.3″ and N 36° 20′ 06.4″, E 050° 58′ 30.7″) from 2,400 to 4,200 m (Figure 1). The two elevational gradients had different extents because of differences in the heights of the peaks.
A total of 133 nested plots were sampled along the two transects (see Table S1). Plots were placed in areas not affected by grazing or other forms of disturbances. In general, there is virtually no human disturbance in the study area, because places are difficult to access and there are few people who visit them for recreation (e.g. mountain climbing). Grazing is very limited and uniform through the whole gradient. We divided each transect into 100 m elevational intervals. In each interval, we took at random three square plots of 100 m2 (i.e. 10 m × 10 m), with the exception of the peak of Mt. Rostam‐Nisht, where only one large plot was taken because of area constraints. From one corner of each 10 m × 10 m square plot, a series of nested subplots (Type I curves) (Scheiner, 2003) of 0.001, 0.01, 0.1, 1 and 10 m2 were sampled for the presence of vascular plant species. We measured the abundance of the plant species per 100‐m2 plot using the Londo scale (Londo, 1976) with the any part system (Dengler, 2009a). For each 10 m × 10 m plot, we recorded elevation and geographical position (with a Garmin 12 XL GPS device) and estimated the percentage of exposed rock visually. At the centre of the same plots, we collected a soil sample (0–10 cm depth) to measure the percentage of sand (by the hydrometer method; Bouyoucus, 1951) and the percentage of total nitrogen (N) as a proxy for soil nutrients (with a Kjeltec System Instrument; Anonymous, 1990). Minimum temperature of the coldest month (MTCM, in °C) and annual precipitation (AP, in mm) were extracted from the CHELSA climate dataset (Karger et al., 2017). CHELSA provides high resolution (30 arc sec) monthly values of temperatures and precipitation for regions with a low density of meteorological stations by downscaling worldwide climate data.
Vector ruggedness measure (VRM, Sappington, Longshore, & Thompson, 2007) was calculated with a 3 × 3 pixel neighbourhood using a digital elevation model with 30 m spatial resolution (AW3D30) provided by the Japanese Aerospace Exploration Agency (JAXA, 2015). Calculations were done with SAGA‐GIS (Conrad et al., 2015). VRM measures the ruggedness of the terrain which is a proxy for habitat heterogeneity. VRM values range from 0 (=flat) to 1 (=very rugged), the latter occurring mainly on ridges or on deep valleys. Finally, we used the Normalized Difference Vegetation Index (NDVI, Tucker, 1979) as a measure of productivity for each vegetation plot. NDVI reflects the photosynthetic activity of vegetation and is therefore particularly useful for evaluating vegetation productivity, health and structure (Pettorelli et al., 2011; Wang, Rich, Price, & Kettle, 2004). For NDVI calculation, satellite imagery of the Landsat 8 mission covering the study area (Path: 165 Row: 35) and sampling period was downloaded from the USGS Data centre as a Level‐1 product. The Level‐1 data were processed using the ‘landsat’ package (Goslee, 2011) and SAGA‐GIS (Conrad et al., 2015). Following the workflow proposed by Young et al. (2017) we converted the raw images to Top‐of‐Atmosphere (TOA) Reflectance using SAGA‐GIS in one step. We did not apply an atmospheric correction as we were not interested in comparing NDVI values from different dates. However, given the mountainous terrain we applied the Gamma topographical correction (Richter, Kellenberger, & Kaufmann, 2009) with a SRTM‐90 m digital elevation model provided by CGIAR (Jarvis, Reuter, Nelson, & Guevara, 2008) to correct the effect of slope and aspect on the reflectance values. All environmental data are given in Table S2.
2.3 Data analysis
We were interested in testing whether the power function of the species–area relationship (SAR) would be the single best model across all elevations or other functions would provide equal or better fits. For this purpose, we tested 18 SAR models for each series of subplots (from 0.001 to 100 m2) in each elevational interval. We used the R package ‘sars’ (Matthews, Triantis, Whittaker, & Guilhaumon, 2019) to fit the various models and to identify the best fit model for each plot. Models which require four parameters, that is, Weibull4 and Beta‐P, were excluded as the AICc cannot be calculated given the small number of samples. The presence of subplots with zero species made the use of the linearized power function problematic because of the logarithmic transformation. Thus, we applied a log (S+1) transformation to all plots. For consistence, we also used S+1 in curvilinear models calculated using the package ‘sars’, but repeated the analyses also with untransformed values of species richness and the nonlinear power model. Models with residuals that deviated significantly from normality and/or homoscedasticity were excluded. We used the log(c)‐ and z‐values obtained from the linearized version in subsequent analyses. Although in a limited number of cases other models gave slightly better fits, we focused on this model in all cases because it provided equivalent parameters to compare (Matthews, Triantis, et al., 2016).
We investigated the relationship between elevation and the power function parameters c and z using a robust regression approach with the ‘rlm’ function in the ‘MASS’ R package (Ripley et al., 2019). Robust regressions are alternative methods to ordinary least squares regressions when data are contaminated with outliers. Robust regressions reduce the impact of outliers by iteratively reweighting least squares. We used a robust regression approach with the ‘MM’ method, which allows finding the optimal weights with a high breakdown point (the proportion of outliers that can be addressed before these observations affect the model). This approach combines the advantage of the M‐ and S‐estimations. The M‐estimation is ‘maximum likelihood type’, which is robust to outliers in the response variable but it is not resistant to the outliers in the explanatory variables (leverage points). The S‐estimation minimizes a robust estimate of the scale of the residuals, and is therefore highly resistant to leverage points. It is also robust to outliers in the response variable, but can be inefficient. The MM‐estimation combines the robustness and resistance of the S‐estimation with the efficiency of the M‐estimation. To adjust the model according to the weights, we used 400 iterations (‘maxit = 400’) with the ‘psi.bisquare’ function for the psi function. We compared the sets of ‘rlm’ models with their AIC values using the ‘MuMIn’ package (Barton, 2017). We focused on the models with delta AIC values ≤2. The same approach was used to test the relationships between the power function parameters and the environmental variables. We conducted separate analyses for elevation and other environmental variables for two reasons. First, while other variables were measured at the level of plot and vary at a local scale, elevation is a variable acting at a broader (landscape) scale. Second, communities are not influenced by elevation per se, but by the several abiotic factors (such as climate, land use, soil, geological settings, etc.) that vary with elevation at the landscape level. Distribution of residuals from rlm‐models was checked using residual versus fitted values plots, QQ‐plots, scale‐location plots and Cook's distance plots. The lack of any trends indicated that the models performed well. To avoid collinearity problems, we first analysed correlations between variables (Table S3). We found a strong correlation of elevation with MTCM and AP in both transects. Rock cover showed a negative correlation with MTCM in both transects, no correlation with AP in transect R and a relatively high correlation in transect L respectively. NDVI was strongly negatively correlated with elevation and positively with MTCM in both transects; negatively with rock cover in transect R and with AP in transect L. Thus, we have decided to not consider AP in the subsequent analyses. The minimum number of parameters was set at 0 (i.e. a model containing only the intercept) and the maximum to 6 (the subsets of possible models with all environmental variables). As the ‘rlm’ function does not provide p‐values, we calculated the significance of each parameter using a Wald test with the R package ‘sfsmisc’ (Maechler, 2019). All analyses were performed in the free statistical software environment R 3.6.1 (R Development Core Team, 2019).
3 RESULTS
Species richness of vascular plants per plot ranged from 0 species (in some 0.001 m2 plots) to 50 species (in some 100 m2 plots) in transect R, and from 0 (in some 0.001 m2 plots) to 43 species (in some 100 m2 plots) in transect L (Table S1). We found large variation in z‐values among the elevational intervals. Namely, z‐values from the log–log model ranged from 0.08 to 0.34 (mean ± SE: 0.22 ± 0.01) in transect R, and from 0.07 to 0.35 (mean ± SE: 0.20 ± 0.01) in transect L (Table S4). Large variations were also found in c‐values from the log–log model, which ranged from 0.15 to 15.93 (mean ± SE: 7.92 ± 0.42) in transect R, and from 1.70 to 16.22 (mean ± SE: 7.22 ± 0.44) in transect L, respectively (Table S4). Both z‐values and log(c)‐values followed the same distribution in the two transects (Kolmogorov–Smirnov tests: D = 0.237, p = 0.052 for z‐values, and D = 0.145, p = 0.502 for log(c)‐values respectively). Values of c and z obtained using the curvilinear model (calculated both with untransformed 0 species values and S+1 transformed values) are also given in Table S4.
The most frequently observed best fit model in both transects was the power function (the best function in 69.74% of cases in transect R, and 66.67% in transect L respectively) followed by the Kobayashi function (18.42% in transect R, and 15.79% in transect L respectively). For the plots in which the power function provided the best fit, the R2 varied from 0.910 to 0.998 in transect R, and from 0.904 to 0.999 in transect L, showing a declining pattern with increasing elevation (Spearman rank correlation coefficient: rs = −0.638, p < 0.0001 in transect R, rs = −0.611, p < 0.0001 in transect L). For the plots in which the Kobayashi function gave the best fit, the R2 varied from 0.771 to 0.993 in transect R, and from 0.914 to 0.996 in transect L, respectively, again declining with increasing elevation in transect L (Spearman rank correlation coefficient: rs = −0.378, p = 0.183 in transect R, rs = −0.826, p = 0.008 in transect L). We found no systematic variation in SAR shape with elevation. In general, plots best modelled by the Kobayashi function were placed at the same elevations where other plots were best modelled by the power function, which typically provided the best fit to two of the three plots, and the Kobayashi function the remaining one. As regards the other models, in transect R the Monod function provided the best fit in only one case, at 3,800 m, whereas the linear model gave the best fit in four cases (at 4,000 m, 4,100 m, 4,300 m and 4,400 m respectively) followed by the power Rosenzweig function (at 4,200 m and 4,300 m respectively) and logarithmic model (at 3,900 m and 4,100 m respectively) in other four cases. In transect L, the linear model provided the best fit in seven cases (at 2,900 m, 3,000 m, 3,400 m, 3,800 m, 3,900 m and 4,000 m respectively) and the Monod function in three cases (at 3,400 m, 4,100 m and 4,200 m respectively) (Table 1). Results obtained using the curvilinear model with 0 species are given in Table S4.
| Mt. Rostam‐Nisht | Mt. Lashgarak | ||||
|---|---|---|---|---|---|
| Elevation | Best model | R 2 | Elevation | Best model | R 2 |
| 2,000 | power | 0.995 | — | — | — |
| 2,000 | power | 0.993 | — | — | — |
| 2,000 | power | 0.996 | — | — | — |
| 2,100 | power | 0.998 | — | — | — |
| 2,100 | power | 0.991 | — | — | — |
| 2,100 | power | 0.998 | — | — | — |
| 2,200 | power | 0.996 | — | — | — |
| 2,200 | Kobayashi | 0.991 | — | — | — |
| 2,200 | power | 0.996 | — | — | — |
| 2,300 | power | 0.994 | — | — | — |
| 2,300 | power | 0.993 | — | — | — |
| 2,300 | power | 0.991 | — | — | — |
| 2,400 | power | 0.995 | 2,400 | power | 0.994 |
| 2,400 | power | 0.998 | 2,400 | power | 0.972 |
| 2,400 | power | 0.982 | 2,400 | power | 0.981 |
| 2,500 | power | 0.987 | 2,500 | power | 0.998 |
| 2,500 | power | 0.997 | 2,500 | power | 0.996 |
| 2,500 | power | 0.997 | 2,500 | power | 0.974 |
| 2,600 | power | 0.985 | 2,600 | Kobayashi | 0.996 |
| 2,600 | power | 0.985 | 2,600 | power | 0.989 |
| 2,600 | power | 0.986 | 2,600 | Kobayashi | 0.993 |
| 2,700 | power | 0.990 | 2,700 | power | 0.997 |
| 2,700 | power | 0.979 | 2,700 | power | 0.986 |
| 2,700 | power | 0.981 | 2,700 | power | 0.988 |
| 2,800 | Kobayashi | 0.991 | 2,800 | power | 0.99 |
| 2,800 | power | 0.982 | 2,800 | power | 0.994 |
| 2,800 | power | 0.989 | 2,800 | power | 0.984 |
| 2,900 | Kobayashi | 0.946 | 2,900 | power | 0.989 |
| 2,900 | power | 0.997 | 2,900 | Kobayashi | 0.990 |
| 2,900 | power | 0.994 | 2,900 | linear | 0.955 |
| 3,000 | power | 0.990 | 3,000 | power | 0.990 |
| 3,000 | power | 0.995 | 3,000 | linear | 0.962 |
| 3,000 | Kobayashi | 0.987 | 3,000 | power | 0.987 |
| 3,100 | Kobayashi | 0.979 | 3,100 | power | 0.983 |
| 3,100 | Kobayashi | 0.993 | 3,100 | power | 0.974 |
| 3,100 | power | 0.974 | 3,100 | power | 0.990 |
| 3,200 | power | 0.984 | 3,200 | power | 0.991 |
| 3,200 | power | 0.995 | 3,200 | power | 0.999 |
| 3,200 | Kobayashi | 0.966 | 3,200 | power | 0.997 |
| 3,300 | power | 0.989 | 3,300 | Kobayashi | 0.993 |
| 3,300 | power | 0.940 | 3,300 | power | 0.981 |
| 3,300 | power | 0.981 | 3,300 | power | 0.989 |
| 3,400 | power | 0.987 | 3,400 | Monod | 0.966 |
| 3,400 | Kobayashi | 0.988 | 3,400 | linear | 0.915 |
| 3,400 | power | 0.993 | 3,400 | power | 0.992 |
| 3,500 | Kobayashi | 0.991 | 3,500 | power | 0.984 |
| 3,500 | Kobayashi | 0.976 | 3,500 | power | 0.974 |
| 3,500 | power | 0.987 | 3,500 | Kobayashi | 0.959 |
| 3,600 | power | 0.998 | 3,600 | Kobayashi | 0.989 |
| 3,600 | Kobayashi | 0.991 | 3,600 | power | 0.940 |
| 3,600 | power | 0.985 | 3,600 | power | 0.905 |
| 3,700 | power | 0.971 | 3,700 | power | 0.976 |
| 3,700 | power | 0.988 | 3,700 | power | 0.967 |
| 3,700 | power | 0.989 | 3,700 | power | 0.940 |
| 3,800 | power | 0.988 | 3,800 | Kobayashi | 0.918 |
| 3,800 | Monod | 0.967 | 3,800 | power | 0.941 |
| 3,800 | Kobayashi | 0.98 | 3,800 | linear | 0.923 |
| 3,900 | power | 0.981 | 3,900 | linear | 0.992 |
| 3,900 | logarithmic | 0.943 | 3,900 | power | 0.904 |
| 3,900 | power | 0.991 | 3,900 | power | 0.922 |
| 4,000 | linear | 0.960 | 4,000 | linear | 0.989 |
| 4,000 | power | 0.97 | 4,000 | linear | 0.976 |
| 4,000 | Kobayashi | 0.963 | 4,000 | power | 0.965 |
| 4,100 | linear | 0.998 | 4,100 | Monod | 0.907 |
| 4,100 | power | 0.968 | 4,100 | power | 0.932 |
| 4,100 | logarithmic | 0.953 | 4,100 | Kobayashi | 0.914 |
| 4,200 | power | 0.932 | 4,200 | Kobayashi | 0.968 |
| 4,200 | powerR | 1 | 4,200 | power | 0.963 |
| 4,200 | power | 0.910 | 4,200 | Monod | 0.957 |
| 4,300 | linear | 0.979 | — | — | — |
| 4,300 | power | 0.978 | — | — | — |
| 4,300 | powerR | 1 | — | — | — |
| 4,400 | Kobayashi | 0.771 | — | — | — |
| 4,400 | power | 0.973 | — | — | — |
| 4,400 | linear | 0.995 | — | — | — |
| 4,500 | power | 0.974 | — | — | — |
We found that z‐ and c‐values decreased significantly with elevation in both transects (Tables 2 and S5, Figure 2). This indicates that, in general, both velocity in species accumulation with sampled area (z) and species richness per unit area (c) decrease with elevation.
| Value | Std. Error | t‐value | p‐value | |
|---|---|---|---|---|
| z | ||||
| R | ||||
| Intercept | 0.416 | 0.026 | 16.107 | |
| Coefficient | −1.000 × 10–4 | 0.000 | −7.694 | 4.583 × 10–11 |
| L | ||||
| Intercept | 0.476 | 0.041 | 11.671 | |
| Coefficient | −1.000 × 10–4 | 0.000 | −6.934 | 7.427 × 10–9 |
| log(c) | ||||
| R | ||||
| Intercept | 3.965 | 0.195 | 20.328 | |
| Coefficient | −6.000 × 10–4 | 1.000 × 10–4 | −10.698 | < 2.2 × 10–16 |
| L | ||||
| Intercept | 4.365 | 0.251 | 17.39 | |
| Coefficient | −8.000 × 10–4 | 1.000 × 10–4 | −10.108 | 3.767 × 10–14 |
In general, z‐ and c‐values were related to the percentage of rock cover, minimum temperature of the coldest month (MTCM) and percentage of total nitrogen (N) (Table 3). Namely, z‐values were positively influenced by MTCM in both transects and by rock cover (albeit not significantly) in transect L (Table 3, Figure 3, see Figures S1–S2 for non‐significant relationships). In both transects, c‐values were positively influenced by MTCM and N, and negatively by rock (Table 3, Figure 3, see Figures S1–S2 for non‐significant relationships). Use of c‐ and z‐values from curvilinear fits gave similar outcomes (Tables S6–S7).
| Intercept | MTCM | N | NDVI | Rock | Sand | VRM | df | logLik | delta | weight | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| z | |||||||||||
| R | 0.377 | 0.011 | 3 | 121.605 | 0 | 0.174 | |||||
| 0.362 | 0.011 | 0.051 | 4 | 122.227 | 0.99 | 0.106 | |||||
| 0.372 | 0.010 | 0.000 | 4 | 121.855 | 1.73 | 0.073 | |||||
| 0.347 | 0.009 | 0.028 | 4 | 121.823 | 1.79 | 0.071 | |||||
| 0.360 | 0.011 | 0.020 | 4 | 121.814 | 1.81 | 0.07 | |||||
| L | 0.435 | 0.017 | 0.000 | 4 | 94.727 | 0.000 | 0.093 | ||||
| 0.364 | 0.016 | 0.000 | 0.001 | 5 | 95.888 | 0.090 | 0.089 | ||||
| 0.348 | 0.015 | 0.001 | 4 | 94.610 | 0.240 | 0.083 | |||||
| 0.353 | 0.013 | 0.068 | 0.000 | 5 | 95.739 | 0.380 | 0.077 | ||||
| 0.426 | 0.016 | 3 | 93.317 | 0.510 | 0.072 | ||||||
| 0.305 | 0.013 | 0.058 | 0.000 | 0.001 | 6 | 96.548 | 1.270 | 0.049 | |||
| 0.353 | 0.015 | −0.022 | 0.001 | 5 | 95.142 | 1.580 | 0.042 | ||||
| 0.370 | 0.013 | 0.046 | 4 | 93.779 | 1.900 | 0.036 | |||||
| c | |||||||||||
| R | 3.278 | 0.092 | 0.665 | −0.005 | 5 | −26.988 | 0.000 | 0.334 | |||
| 3.442 | 0.091 | −0.006 | 4 | −29.129 | 1.990 | 0.124 | |||||
| L | 3.707 | 0.131 | 0.381 | −0.005 | 5 | 4.001 | 0 | 0.284 | |||
| 3.991 | 0.136 | 0.409 | −0.005 | −0.004 | 6 | 4.747 | 1.010 | 0.171 | |||
| 3.412 | 0.114 | 0.342 | 0.261 | −0.005 | 6 | 4.451 | 1.60 | 0.127 | |||
- Abbreviations: Delta, difference from each model to the best model; df, degrees of freedom; logLik, Log likelihood value; MTCM, minimum temperature of the coldest month (°C); N, percentage of total nitrogen (%); NDVI, Normalized Difference Vegetation Index; Rock, percentage of exposed rock (%); Sand, percentage of sand (%); VRM, vector ruggedness measure; Weight, Akaike weight showing support of each model.
4 DISCUSSION
There is a long debate about which model should best describe the SAR and numerous functions have been proposed (Dengler, 2009b; Tjørve, 2009; Triantis et al., 2012). Our multi‐model analysis revealed that the power function, followed by the Kobayashi function, provided the best fit to most of our nested plots. The power model has been reported as the best fit model in many comparative studies (Arrhenius, 1921; Azovsky, 2011; Dengler, 2009b; Dengler & Boch, 2008; Drakare et al., 2006; Matthews, Guilhaumon, et al., 2016; Rosenzweig, 1995; Triantis et al., 2012; Williamson, 1988), followed by the Kobayashi function (Triantis et al., 2012). In general, for the vast majority of the elevational belts, the best‐performing models have two parameters with a convex shape (Matthews, Guilhaumon, et al., 2016; Matthews, Triantis, et al., 2019; Tjørve, 2009; Triantis et al., 2012). Exceptions are the linear function, which has a linear shape, and the power Rosenzweig function, which has three parameters with a convex shape. The biological meaning of the parameters of most functions proposed to model the SAR is, however, obscure or controversial (Dengler, 2009b; Triantis et al., 2012). In general, it is expected that the models with more parameters will be more flexible to fit the datasets, but our results indicated that simpler models provided the lowest AIC values in most cases, as already reported in other studies (e.g. Triantis et al., 2012), which indicates that the extra flexibility of more complex models is not ‘worth’ the cost of their extra parameters. Also, Tjørve (2009) expressed the idea that models with fewer parameters, such as the power and Kobayashi functions, are often more preferable when there are few data points and large scattering, whereas more complex and flexible models, such as the power Rosenzweig function (which is a power function modified with the addition of a constant to improve the fit), may give better results when simpler models do not fit well.
An advantage of the power function is that it is controlled by two parameters which have been largely investigated and can be easily interpreted in an ecological context (Fattorini et al., 2017; Matthews, Guilhaumon, et al., 2016). The parameter c represents the expected mean number of species per unit area, whereas z can be interpreted as a scaling factor describing how fast the response of species richness to area changes along the SAR curve. As a general rule, it has been suggested that z should increase with isolation, species trophic ranks, nestedness and spatial aggregation of the individuals, and should decrease with species dispersal ability, abundance of common species, human impact, latitude (possibly as a response to increasing energy availability) and elevation (Fattorini et al., 2017; Qiao, Tang, Shen, & Fang, 2012). The z‐parameter also varies with the scale. For example, Fridley, Pett, Wentworth, and White (2005) found that z was high at fine scales, lower at intermediate scales and again rising at large scales. Triantis et al. (2012) also found that z‐values appeared to vary in relation to the range of island areas encompassed: for island datasets spanning just two orders of magnitude the mean value of z was significantly higher than for datasets spanning more orders of magnitude of island area. Simulations indicate that z‐values increase with increasing extinction rates and decreasing colonization rates (Rybicki & Hanski, 2013). However, landscape disturbance and fragmentation may increase z even without species loss, if the spatial distributions of species become more localized after disturbance, or the relative species abundance distribution (SAD) becomes more skewed (He & Hubbell, 2013).
Elevational patterns in c‐values found in our study showed that the highest values of log(c) were at the lowest elevations, while the highest elevations encompassed the lowest values. This is an indication that species richness per unit area decreased along the elevational gradient (see also Moradi & Oldeland, 2019). A continuous decline of species richness with increasing elevation is considered the most widespread pattern of diversity variation along elevational gradients, although there is substantial evidence for the existence of mid‐elevation peaks in a broad range of organisms (Lomolino et al., 2010; McCain, 2010; Rahbek, 1995; Sanders & Rahbek, 2012; Stevens, 1992). However, our study involved a section of the gradient that starts well beyond the mid elevations. Overall, the ecological space involved in our study encompasses variations in abiotic factors (such as the climatic capacity for productivity), and biotic conditions (such as the species pool and species interactions), which can offer the most obvious explanations for the monotonic decrease in species richness with elevation (Kaspari, O’Donnell, & Kercher, 2000; Tjørve, 2009). However, the commonly observed decrease in species richness at higher elevations might be also a possible reflection of a reduction in available area because of the conical shape of mountains (McCain, 2010). The use of c‐values of SARs constructed at different elevations can overcome this problem. Since c‐values express the expected mean number of species per unit area, they represent estimates of species richness not biased by different extent of the mountain surface at different elevations. It is important to note that using c‐values from the SAR is a more correct approach than calculating species density measures (number of species per plot area). Species richness does not increase linearly with area, thus it is not correct to calculate species richness per unit area by simply dividing the recorded number of species by the size of the sampled (plot) area. Since it is always subjective to choose the size of a sampling plot to record species richness, the use of plots of different size would produce results that are not comparable. Thus, taking into account that richness does not vary linearly with area, the use of a SAR parameter that expresses the expected mean number of species per unit area across variation in sampling areas makes this type of measure of richness more straightforward than measures of species density calculated for any given predefined sampling area. On one hand, differences in c‐values can be interpreted as a direct result of differences in the richness of the regional species pools. Since the local species richness can be perceived as samples drawn from a regional species pool, the species richness recorded at the scale of the plot can be influenced by the overall number of species present at the elevation from which a plot has been censused (Austrheim, 2002), which, in turn, can be a reflection of the area available at that elevation (Romdal & Grytnes, 2007). Therefore, the low c‐values recorded at high elevations might be a consequence of the smaller number of species present in the higher vegetational belts, as a result of the smaller available surface and less favourable environmental conditions. In the Alborz Mountains, higher elevations with harsher environments filter out species by imposing challenging conditions to plant life, effectively reducing the species pool in the area (Moradi & Oldeland, 2019). Yet, it is not difficult to imagine situations where different mean numbers of species per unit area can be found in sets of areas that have globally the same overall species richness, but which differ in SAD and species spatial distributions, two factors that influence the probability of a species of being included in the samples (He & Legendre, 2002; Tjørve & Tjørve, 2008). For example, two belts might have the same overall richness, but different c‐values, if they differ in the number of rare species (i.e. species with small relative cover), because these have a lower probability of being included in the plots, and thus will be under‐represented, leading to smaller values of species richness and c‐values.
In the power function, z is related with β‐diversity (species turnover) in different ways, according to the type of SAR. In fact, the relationship between z and β‐diversity is not obvious, and z cannot be considered a true and direct measure of β‐diversity, as defined by Tuomisto (2010), although it indicates differences in species composition. However, in SARs arising from nested plots (Type I in Scheiner's, 2003 classification), z reflects directly species gains with increasing surface, whereas other types of SARs (Type II, III and IV in Scheiner's classification) involve different forms of averaging of β‐diversity, as defined by Scheiner (2003).
We found that z‐values decreased with elevation, which indicates a decrease in species accumulation, possibly a reflection of the smaller ‘regional’ pool (i.e. a smaller number of species present in the highest elevational belts).
Despite many studies on the SAR, factors and mechanisms underlying variation in the c and z parameters are still unclear (Triantis et al., 2012). It is difficult to evaluate the relative importance of different environmental factors in determining how these parameters vary along elevation, since several environmental variables change directly and indirectly with elevation (Körner, 2007). Climatic variables were strongly correlated with elevation in our study, yet rock cover explained an important fraction of variation in z‐ and c‐values which is not explained by climate or elevation alone. Climatic stress was already known to strongly influence plant diversity along elevation in the Alborz Mountains (Moradi & Oldeland, 2019). In the present study, we found that temperature (MTCM) explained a large amount of variation in c‐values. Namely, the positive relationship between MTCM and c‐values, and the inverse correlation between MTCM and elevation, indicate that colder climatic conditions at higher elevations have a negative impact on plant diversity. However, as elevation increases towards the nival belt, not only the associated climatic severity changes but also resource restrictions limit plant growth (Körner, 2007). Soil characteristics, including nutrients and soil structure, are known to be among the variables affecting diversity, albeit with less importance when compared to other physical variables (Austrheim, 2002). In our study, at high elevations, in addition to climatic stress, the increase in rocky surface limits plant cover to scree vegetation (Moradi & Oldeland, 2019), which reflects a certain type of habitat where plant colonization is limited and only specialists can grow. Rock cover reduces the available niche space and increases the pressure on species to adapt to this challenging environment (Moradi, Attar, & Oldeland, 2017). As a result, high rock cover causes a less dense vegetation at high elevations, that is, fewer individuals, and hence fewer species per unit area. Nitrogen had a positive influence on specie richness, which suggests that a higher nutrient availability promotes diversity. Temperature had a positive influence also on z‐values, indicating that more favourable conditions promote a higher species accumulation rate, thus paralleling the results obtained for the c‐values.
Previous research found a negative relationship between productivity (measured as biomass) and z‐values in plant plots, whereas the relationship with c‐values was found to be either positive or absent (Chiarucci, Viciani, Winter, & Diekmann, 2006; Pastor, Downing, & Erickson, 1996). We found no strong influence of NDVI on either c‐values or z‐values, although weak positive trends have been found in both cases. Positive relationships of both c and z with productivity point to simultaneously high local and high spatially distributed richness in productive systems, that is, productivity and diversity will be simultaneously maximized (Pastor et al., 1996). These results, coupled with the decreasing behaviour of NDVI, c and z with elevation, suggest that high altitude vegetation is less diversified and less productive.
5 CONCLUSIONS
Exploring variation in SAR shape along an elevational gradient allowed us to investigate the meaning of the power function parameters from a new perspective. We found that the c‐values decreased with elevation, which is consistent with the commonly reported altitudinal decrease in species richness. This reduction in species density is explained by the increase in bare rock. Minimum temperature and percentage of total nitrogen affected positively c, thus indicating that these two factors that promote plant growth also promote plant diversity. Similarly, z‐values decreased with elevation, and were positively influenced by temperature, which indicates that progressively harsher conditions decrease species accumulation rates. However, our results are based on two transects of the same mountain system and the relative importance of environmental characteristics can change according to the overall geographical setting, which calls for further research in other mountain systems.
ACKNOWLEDGEMENTS
We express our thanks to the Vandarbon Mountaineering Federation staff and to the inhabitants of Delir rural district for their assistance during fieldwork. The generosity and cooperation of local shepherds made this study possible. We thank them for their interest in our study and for their hospitality and kindness. We are grateful to Ole Vetaas, Tom Matthews and an anonymous reviewer for their very constructive comments on a previous version of the manuscript. We thank the DAAD (Deutscher Akademischer Austausch Dienst) for a grant to H.M. (funding grant number 57214227). No permits were required for field work.
CONFLICT OF INTEREST
We declare no competing of interests.
REFERENCES
BIOSKETCHES
H. Moradi is interested in vegetation ecology, plant functional and taxonomic diversity and patterns of diversity. She studies both vascular plants and bryophytes with a focus on mountain forests and alpine grasslands.
S. Fattorini has broad interests in island biogeography, conservation biogeography, macroecology, community ecology and urban ecology. He is especially interested in coastal, desert and montane environments.
J. Oldeland is a vegetation ecologist with a strong background in ecological remote sensing and biostatistics. His further interests are biodiversity, drylands and high mountain ecosystems.
Author contributions: H.M., S.F. and J.O. conceived the ideas; H.M. collected the data; H.M., J.O. and S.F. analysed the data; H.M. led the writing; S.F. and J.O. contributed in the form of discussions and suggestions. All authors approved the final manuscript.
